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■ Finding the distribution of vibro-acoustic energy in complex built-up structures in the mid-to-high frequency regime is a 
CN • 
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difficult task. In particular, structures with large variation of local wavelengths and/or characteristic scales pose a challenge 
referred to as the mid-frequency problem. Standard numerical methods such as the finite element method (FEM) scale with 
the local wavelength and quickly become too large even for modern computer architectures. High frequency techniques, such 
as statistical energy analysis (SEA), often miss important information such as dominant resonance behaviour due to stiff 
' or small scale parts of the structure. Hybrid methods circumvent this problem by coupling FEM/BEM and SEA models in 
a given built-up structure. In the approach adopted here, the whole system is split into a number of subsystems which are 
treated by either FEM or SEA depending on the local wavelength. Subsystems with relative long wavelengths are modeled 
using FEM. Making a diffuse field assumption for the wave fields in the short wave length components, the coupling between 
subsystems can be reduced to a weighted random field correlation function. The approach presented results in an SEA-like 
set of linear equations which can be solved for the mean energies in the short wavelength subsystems. 
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> ' 1. Introduction 



Determining the distribution of vibro-acoustic energy in complex built-up structures using finite element methods (FEM) is 



in 

, often a difficult task, in particular in the high frequency limit. The number of degrees of freedom (DoF) in FE models scale 
with the wavelength of the underlying wave field and computations quickly become prohibitive with increasing frequency. 
Moreover, "numerically exact" results obtained by an FE approach for a specific structure may be of little practical value. 
' The vibro-acoustic response of "identical" structures assembled, for example, as part of a manufacturing process, is very 
sensitive to small changes in material parameters and/or variability of the shape of the structure. These differences may lead 
to large changes in the resonance spectrum and most of the detailed information of a full (and expensive) FEM calculation 
for an individual example has then at best statistical significance. In room acoustics, this computational ceiling may be as 
low as a few 100 Hz for moderate room sizes, it is roughly at 200 Hz for FE models describing the interior noise levels in 
^ ! cars [HE]. 

Statistical Energy Analysis (SEA) [3] is a popular tool to tackle the mid-to-high frequency regime. In an SEA framework, 
the whole system is viewed as a number of interacting subsystems giving rise to a set of linear equation for the mean wave 
energies stored in each subsystem. The SEA approach assumes that the wave field in each component is diffuse thus having 
universal statistical properties (that is, properties not depending on the shape, material parameters or even the wave equation 
characterizing a specific subsystem). SEA results may be interpreted as an average response of an "ensemble" of similar 
structures or of an individual system after averaging over frequency intervals large on the scale of the mean resonance density. 
SEA modelling is particularly useful if the wavelength is much smaller than characteristic scales of individual subsystems 
across the whole structure. 

In practical applications, this assumption often does not hold. Local wavelengths may vary considerably across a structure 
and may at places become comparable to the size of a subsystem. One may, for example, encounter stiff or narrow parts 
having only a few resonances in a given frequency band which are coupled to subsystems with a high modal density where SEA 
assumptions are valid. The long wavelength components will give rise to characteristic resonance features which dominate 
the vibro-acoustic response - but can not be reproduced in an SEA approach. This leads to mid-frequency bands where the 
applicability of FEM and SEA do not overlap. To overcome such problems, a range of hybrid methods has been developed 
based on the concept of distinguishing between deterministic (stiff) and stochastic (short wavelengths) components. 
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Figure 1. (Colour online) The system is split into subsystems: p\ andp2 are stochastic subsystems, ri\ and rii are deterministic 
subsystems. Deterministic boundaries are shown as full lines and random boundaries as dashed lines. Dotted lines denote 
interfaces between subsystems. 

FEM is the method of choice for the deterministic subsystems while for the stochastic subsystems, one makes simplifying 
assumptions which reduce the computational cost and make it possible to find the response averaged over an ensemble of 
statistical subsystems. One such hybrid method is the fuzzy structure theory [H [5] which predicts the response of a master 
structure coupled to a large number of secondary structures. Assuming randomness in the secondary structures, the effects 
of the fuzzy substructures on the master structure can be described by random impedance operators. This approach both 
reduces the number of DoF and gives an average response. Other hybrid methods include, (i) analytical impedance techniques 
[51(7], where floppy parts of the structure are modelled as receiver impedances; (ii) modal |8] and SEA/modal [S] approaches; 
and (iii) the "Smooth Integral" /FEM hybrid method [TU], where probabilistic formulations are introduced as a modification 
of a boundary integral method. 

A hybrid method based on wave concepts matching the solutions in the two kinds of subsystems without a priori 
assumptions about mode coupling was first developed by Shorter and Langley [111 I12j . The central idea of the method 
is a reciprocity result regarding the forces exerted on the boundaries of the deterministic subsystems |13j . It is assumed that 
the wave field in stochastic subsystems can be decomposed into "direct" and "reverberant" field components and that the 
reverberant field component is described in terms of a random diffuse field. The reciprocity relationship couples the direct 
field radiation with diffuse reverberant loading at the interfaces between the deterministic and stochastic subsystems. 

In this paper, we also divide the whole structure into stochastic and deterministic subsystems and split the global wave 
field into a direct and a reverberant component. We then introduce two important modifications from the Shorter and 
Langley approach 112) . Firstly, we compute the direct field explicitly and coherently across the whole structure - details 
will be provided in Sec. [31 This yields a direct coupling between deterministic subsystems which is important in particular 
for medium to strong damping. Secondly, we describe the response of deterministic subsystems to external excitations in 
terms of Green functions. The coupling between the reverberant field in the deterministic and the stochastic subsystems is 
then facilitated by a diffuse field-field correlation function acting on the Green functions of the deterministic subsystems. 
The resulting expressions are formally equivalent to the reciprocity relation derived in |13j using the force-force correlation 
function as has also been pointed out in [14] . However, our approach simplifies the derivation considerably and offers an 
alternative viewpoint regarding possible future advancements of the method. The contributions of both wave fields will be 
finally coupled through energy balance equations. 

The article is organized as follows: the general set-up is introduced in Sec. [5] In Sec. [3J we discuss a method to compute 
the direct field. In Sec. 01 the average response of a deterministic subsystem to an energy influx from adjacent stochastic 
subsystems is described in terms of a random field correlation function. In Sec.[5j we collect results from the previous sections 
and present the outline of the method. Numerical examples are given in Sec. H3 

2. The general set-up 

We will work with the Helmholtz equation in two spatial dimensions here, though the concepts are valid for general wave 
equations such as vectorial elasticity as well as for 3D problems. We will compute the solution to the wave equation for a 
set of coupled acoustic 2D cavities. To stay within an acoustics framework, we will use Neumann boundary conditions on 
the boundaries of the domains, although any other linear boundary condition could be used instead. 
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The hybrid method is now set up as follows: the whole system is divided into deterministic and stochastic subsystems. 
Deterministic subsystems are those that, at a given frequency, are characterized by a wavelength comparable to the size of 
the subsystem. In actual applications, these may be bolts, frames, narrow joints or wave guides to name a few. Stochastic 
subsystems are those for which the local wavelength is small relative to typical scales of the system. The wave field is 
generally characterized by rich interference patterns and its statistics is often well described by that of a Gaussian random 
field or generalizations thereof |15l 116] . We assume that we know the geometry of deterministic subsystems exactly and call 
the corresponding boundaries deterministic boundaries. For stochastic subsystems, we will not specify the actual form of the 
boundaries, only the area of the subsystem will enter the equations. We will often consider ensemble averages over stochastic 
subsystems by changing the boundary (keeping the area fixed) - consequently, we will call the boundaries of stochastic 
subsystems random boundaries, see also [13] . The ensemble average can be considered as averaging over a range of 'similar' 
systems, such as cars produced on an assembly line differing by small changes in fabrication. 

In the following, we assume that deterministic subsystems are connected to stochastic subsystems only. If two determin- 
istic subsystems share a boundary, they are merged into one deterministic subsystem. The same holds for the stochastic 
subsystems. The whole system then consists of N deterministic subsystems and P stochastic subsystems, say. Fig. [1] shows 
a particularly simple example with N = P = 2; this system will be later considered for a numerical validation of the method. 
Our aim is to find the mean wave energies in the stochastic and deterministic subsystems, (E p ), p — 1, . . . ,P and (E n ), 
n = 1,...,N, respectively. 

We are looking here for approximate solutions of the stationary Helmholtz equation satisfying Neumann boundary con- 
ditions on all boundaries, that is, 

o- V 2 ^ + 0> 2 - = f, (1) 

where ip denotes the acoustic pressure field over the whole system, 77 is a damping coefficient, a is the stiffness, p is the 
density and / is a forcing term leading to the excitation of the system. In what fallows we set p = 1. We will work with 
dimensionless units throughout; in the particular example, Fig. [TJ the width of the deterministic subsystem ni,Ti2 is scaled 
to one. 

We now decompose the global wave field into a direct and a reverberant field across the whole structure, that is, we write 

ip = tp W + </> (r) . (2) 

Both these fields are defined everywhere in the structure, that is, both in the deterministic and stochastic subsystems. The 
direct field ip^ satisfies the wave equation ([TJ) with Neumann boundary conditions on deterministic boundaries and absorbing 
boundary conditions on random boundaries. It thus corresponds to a wave field being reflected from and scattered by the 
boundaries of deterministic subsystems only, (full lines in Fig. [1}. Random boundaries, on the other hand, are ignored, that 
is, the system is considered open at these boundaries for the purpose of calculating ipw. For an example of a direct field, 
see Fig. H 

The stochastic reverberant field tp^' is the part of the total field ip which originates from multiple reflections on both 
the random and deterministic boundaries. It can be considered as the wave field emerging after the first reflection of ip^ 
at the random boundaries. The boundary conditions for ip^ are implicitly defined through the boundary conditions of 
the wave equation and value of the field ip^ on the random boundaries. As we will see later, there is actually no need 
to calculate ip^ explicitly in our hybrid method and results will be independent of the boundary conditions chosen on the 
random boundaries. Ensemble averages will be performed by changing the boundaries of the stochastic systems leaving the 
area invariant; from the set-up, it is evident that these changes will only affect ip^ r \ whereas the direct field ip^ remains 
independent of the shape of the random boundaries. The reverberant field, on the other hand, is extremely sensitive to 
changes of these boundaries. Due to the size of the stochastic subsystem, small deformations of the boundaries will lead to 
large-scale changes in the interference patterns of ip^ originating from the multiple reflection processes. This opens up the 
possibility to describe the reverberant field in a statistical sense. The direct field is thus independent of the reverberant field 
ip( r \ but acts as a source term for the reverberant contribution ip^ r \ 



3. Direct field 

We start by determining the direct field contribution ip^ of the global solution ip in Eq. The areas of the deterministic 
subsystems have been meshed to obtain FEM models for each subsystem. A part of the structure in FigQ] together with 
the mesh is shown in Fig. [2] in more detail. Here, one can see the boundary of the structure, the mesh used in the FEM 
model of the deterministic subsystem and the interface which is introduced to separate the deterministic and the stochastic 
subsystem. We have a certain amount of freedom to chose this interface, that is, changing the form of the interface within 
the stochastic subsystem will not change the overall result. Notice that in Fig. [5] the FEM mesh is somewhat expanded into 
the stochastic subsystem in contrast to the interface shown in Fig. Q] 
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3.1. Dynamic stiffness matrix for a single deterministic subsystem 

To start with, let us assume that there is only a single deterministic subsystem, that is, we neglect any interaction between 
different deterministic subsystems. According to the definition given in the previous section, the direct field then satisfies 
Neumann boundary conditions on the deterministic boundaries and transparent boundary conditions on the interfaces with 
the neighboring stochastic subsystems. Such a model accounts for the direct field radiation into stochastic parts, that is, 
outgoing waves emanating from interfaces of the deterministic subsystem will propagate outwards and will not return to the 
subsystem . Describing the wave problem for the deterministic subsystems in an FEM setting, the wave equation takes the 
form of a set of linear equations 

flWqW =g (n) , (3) 

where n is the deterministic subsystem under consideration. Here, q is the vector of nodal DoF, g is the vector of nodal 
forces, and D is dynamic stiffness matrix set up as usual in terms of stiffness, mass and damping matrices. 

There are various ways to implement transparent (absorbing) boundary conditions on the interfaces. One of the possible 
approaches to achieve absorbing boundary conditions is to construct an operator relating Dirichlet data to Neumann data on 
the interface. Such Dirichlet-to-Neumann map |17j can be cast in the form of what has been recently reported as hybrid FE 
- wave based technique in jTS]. In particular, the application for a two dimensional unbounded problem is explained in Ref. 
[19] . In what follows we will use the method from Ref. [TI5] to construct the dynamic stiffness matrices of the deterministic 
subsystems. However, other methods such as perfectly matched layer (PML) techniques 20, 21 could also be used. 

In the following it will be advantageous to split the equation of motion for the deterministic subsystem into deterministic 
DoF qd and interface DoF q,, that is, 

/ r\M n(") \ / _.(«) \ / _(«) \ 

(4) 

where the subscripts d, i denote DoF which lie inside deterministic subsystems and on interfaces, respectively. In particular, 
the sub-block Di accounts for absorbing boundary conditions for outgoing waves. In the setting considered here, the sources 
are in the deterministic subsystem and the DoF on the interface, q^, only contain outgoing wave components. We now move 
on to the case of several deterministic subsystems coupled through a direct field. 

3.2. The direct field propagator 

In a next step, we will describe the coupling between different deterministic subsystems in terms of free wave propagators 
mapping outgoing waves from one deterministic subsystem into incoming waves on another deterministic subsystem. In the 
following, we will use the notation (np) for the interface between the n t h deterministic subsystem and the pth stochastic 
subsystem assuming n is adjacent to p. We denote by q|™ p ^ the direct field components on the (np) interface; the corresponding 
points are denoted with 'stars' in Fig. O We will distinguish between incoming and outgoing components with respect to 
the deterministic subsystem on the (np) interface and write 

Mi — Hi n T Hout ■ l°J 

We will leave out the superscript (d) (for direct field) in this section - all field quantities considered will be components of 
the direct field. 

We also restrict ourselves to circular interfaces. This choice is not a severe limitation of the method and allows for very 
general configurations due to the freedom we have to choose the interface between adjacent components. Let us now assume 
that we have found the wave solution in the deterministic subsystems. In particular, let us assume we know the outgoing wave 
component on the (np) interface. (The actual coupling of FEM models corresponding to different deterministic subsystems 
will be described in the next section.) The direct field can now be extended into the stochastic subsystems. Due to the 
circular shape of the interface, we expand the outgoing field, qi™^, in terms of a Fourier basis, 

(j>W(6) = sin(2TO(9), (j} [ m ] (e) = cos(2to(9), (6) 

where 9 is the angle parameterising the interface. That is, we write the outgoing wave field as 

1 9 M 

^(0) = -c + -J2 (<4 + W(*) + c$-wvj) (?) 

m=l 

with expansion coefficients 

c (n P ) = [Coi c[ + \c{-\ . . . , cfiKctff. (8) 
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Figure 2. (Colour online) a) An (np) interface: solid lines are the boundary of the system, dashed lines, the interface between 
stochastic and deterministic subsystems. The interface is chosen in form of semicircle centred at rg (marked by a cross), b) 
The meshing for FEM computations: stars denote the degrees of freedom on the interface with field components q^, circles 
denote the degrees of freedom in the vicinity of the boundary inside the deterministic subsystem with field components d. 



The wave field ip^t W corresponds to the wave vector q^f in Eq. ([5]) on the interface mesh points. The coefficients C^ np ^ 
can be obtained through numerical integration. The latter may be written in the form 



c (n P ) = ytaOqfcri 



with VC 



(np) 



;)A ; 



(9) 



The matrix V performs the numerical integration over the interface to obtain the Fourier coefficients of the field q^ on the 
nodes of the FEM mesh. Here, the index I labels the degrees of freedom on the interface and A; is the summation weight of 
the quadrature. The truncation parameter M is determined by the wave length of the problem. 

The direct field emanating from the n t h deterministic subsystem into the pth stochastic subsystem through an (np) 
interface can now be write explicitly as 



V> (np) (r) = 



Co 



H^(kr) 



M 

E 

m— 1 



n (+) A+) , n 



(10) 



where Hm are m t h order Hankel functions outgoing with respect to the deterministic subsystem, r = (r, 9) are the cylindrical 
coordinates centred at tq (cross in Fig. [T] a) and R is the radius of the interface. Note, that the wave number k can be 
complex in case of non-zero damping. The solution written in the form of Eq. (1 1 0[) corresponds to radiation from the 
deterministic subsystem into the stochastic region; it will be used to evaluate solutions only on the interface and outside of 
the nth deterministic subsystem, that is, for r > R. Of course Eq. (JTUJ) implies that the stochastic subsystem is homogeneous 
in a sense that the material parameters are constants. 

The direct field emanating from the deterministic region n into the stochastic subsystem p adjacent to n is fully described 
by the vector C. The field at a point r inside the stochastic subsystem p can now be obtained using Eq. (fT0|) . That is, we 
obtain 



1 /,("P)(r) = L {np \T) 



„(™p) 

Hon* ' 



where 



and 



£( n P)(r) = £("P)( r )y( np ). 



H£l(kR) 



(11) 
(12) 

(13) 
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In what follows, iP 1 "^ is termed the (np) component of the direct field propagator. 

To proceed we introduce a set of direct field propagators which are defined implicitly as 

c&' p) =L^c&*\ (14) 

for n 7^ n' . That is, L tynpn ) denotes the direct field propagator mapping outgoing waves on the (np) interface onto incoming 
waves at the (n'p) interface where the deterministic subsystems n and n' are connected through subsystem p. Note, that 

q[™ p ' corresponds to incoming wave contributions at the (n'p) interface for n ^ n' . We define L^ npn ^ = which means that 
the direct field propagator produces only incoming waves with respect to the deterministic subsystems. We can now define 
the total direct field propagators L^ p ' corresponding to the stochastic subsystem p in the form 



/ ... L^ 1 ) 



V i< lpAr ) • • • o 



(15) 



Note that contains only those direct field propagators ]J^ lpn ' with n and n' adjacent to the stochastic subsystem p. A 
global direct field propagator is now defined as 

1=1 I. (1(i) 




L (P) 

Having constructed the global direct field propagator, Eq. (JT5J), one can related the incoming to the outgoing field components 
on all interfaces, that is, 

q in = Lq out . (17) 
3.3. The global equation of motion for the direct field 

We saw in the previous subsection that the direct field in the stochastic subsystems is determined by the solution on the 
interfaces alone, as we do not account for any reflections from random boundaries. Our treatment thus automatically fulfills 
the free radiation (or absorbing) boundary conditions on the random boundaries. We proceed by constructing the direct field 
solutions in the deterministic subsystems and on the interfaces on an FE mesh making use of the direct field propagator. 
Using Eq. (|17[) . one can couple Equations of motion Q for individual deterministic subsystems to obtain 



D d D dl (I + L) \ / q d \ = / g d 
Did A J V <w / V Si 



(18) 



where / is the identity operator on the boundary DoF. The matrix Dd is a block diagonal matrix set up by the sub-matrices 
D d n ^ defined in Eq. Q and similarly for Ddi, Did, Aj that is, all matrices and vectors in (TT5)) are defined globally across all 
deterministic subsystems and their interfaces. Note, that Di accounts for absorbing boundary conditions for the outgoing 
fields at interfaces. Incoming waves responsible for the coupling between wave fields in different deterministic substructures 
are taken into account through the direct field propagator L via Eq. (jTT]) . 

Equation (1181) is now solved for and (\ out . The incoming component q„ can then be obtained from Eq. (|17|) . The 
solution of Eq. (TT5|) gives the wave field produced inside the deterministic subsystems and on the interfaces including coupling 
of different deterministic subsystems via radiation flowing through stochastic subsystems. Damping in stochastic subsystems 
is taken into account by choosing complex wave numbers in the direct field propagator L. Note that the direct field is 
coherent across the whole structure and that the deterministic subsystems exchange energy through the direct field. If we 
have only one deterministic subsystem as considered in Sec. 13.11 the method introduced is an implementation of absorbing 
boundary conditions on the circular interfaces as for example considered in [19, 22 . The direct field can be found everywhere 
in the stochastic region with the help of the direct field operator L defined in Eqs. (1111) . (fT5"]) . ([To]) . A solution for the direct 
field in the structure shown in Fig. [1] is presented in Fig. [3] It can clearly be seen that the direct field propagates freely 
through the random boundaries. 

We have so far assumed that the structure is excited in a deterministic subsystem. If the forcing is performed inside a 
stochastic subsystem p, we have to add to the right hand side of Eq. (|18p a force acting on the (np) interfaces by the incoming 
wave field emanating from the source. The latter can be calculated using the free field Green function inside p. 
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Figure 3. (Colour online) The logarithmic plot of a solution for the energy in the direct field. 



In order to compute the energy flux transported via the direct field from deterministic to stochastic subsystems one needs 
the direct field components and the total field on the boundary q 4 . They can be obtained from the solution of Eq. (|Tg)) 
together with Eq. (jf 7[) . The direct field energy input into the pth stochastic subsystem per period T — 2*k /lj is then given as 

N 

p l d) = f £ ^r ] n {d^} q^i, (19) 

71=1 

where D^ p ^ is sub-matrix of Did which account for interactions at the (np) interface. In fact, to use Eq. (fT9)l one only needs 
to evaluate qjj at mesh points inside n next to the interface (marked by circles in Fig. [2]). 

An extension of the method to 3D can be performed along the same line. In this case, the interfaces have the form of 
spherical segments. The segments of a sphere defining the interface can be of fairly arbitrary shape. Finding a 2D, orthogonal 
basis on these interfaces in analogy to the Fourier basis in Eq. © will lead to an additional numerical overhead. However, 
the general principle remains the same. 

4. The stochastic reverberant field 
4.1. FEM equations 

To account for the response of the stochastic subsystems, we need to estimate the stochastic reverberant field ip^ in addition 
to the solution of Eq. (|18[). The stochastic field is created in interference of waves reflected from the random boundaries. In 
what follows, we will assume that ipW is a random diffuse field, see |15] and references therein. We want to avoid an exact 
computation of 4> either because we do not have enough information about the random boundaries or because a numerically 
exact treatment of the stochastic systems would need exorbitant large FEM models. In our hybrid approach, we will thus 
only take into account the response of a deterministic subsystem to the reverberant field exited in the neighboring stochastic 
subsystems. We will furthermore assume that the reverberant fields in different stochastic subsystems are uncorrelated and 
that the reverberant field inside the n t h deterministic subsystems is an incoherent superposition of the partial fields excited 
through the different interfaces of n. It is thus sufficient to compute the response of the nth deterministic subsystem to a 
diffuse excitation at every (np) interface separately for all p adjacent to n. 

The stochastic field in subsystem p leads to stochastic field components on the (np) interfaces which can be modelled in 
terms of stochastic reverberant forces. These forces in turn excite vibrations in the deterministic subsystems n and lead to 
transport of vibro-acoustic energy through subsystem n to other stochastic subsystems p' adjacent to n. We will denote the 

reverberant field in the n t h deterministic subsystem as [q^,q- r ^] and the discussion will be restricted to the local field in 
n from now on. We will also omit the superscript (r) for reverberant field in this section - all field quantities considered will 
be reverberant. 
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We start by computing the reverberant field in a deterministic subsystem n excited on the (np) interface by a stochastic 
field in subsystem p. We again split the solution on the interfaces (here of subsystem n only) as 



(n) (n) (n) 



and write the FEM equations in the form 




f(np) 





(20) 



with the same dynamic stiffness matrices as defined in Eq. (jj). It implements absorbing boundary conditions at the interfaces, 
that is, the interfaces act as openings through which outgoing waves leave the system without being reflected. The right 
hand side of Eq. (|20p contains the stochastic reverberant force 

f (np) = _£><>P) q W (21) 

acting on the deterministic DoF next to the (np) interface (marked by circles in Fig. [5J . Note that the only non-zero entries 
of the force vector f( n w are those along the (np) interface. We of course do not know qi„, so this variable needs to be 
eliminated eventually. 

The solution of Eq. (|2U)) for a given force term can be found with the help of the Green function G^™ p ' acting on the DoF 
inside subsystem n with source points on the (np) boundary. It is the solution of the equation 




A n p) 



(22) 



where Jo ^ is the identity matrix on the DoF next to the (np) interface marked by circles in Fig. [2] 

The Green function in Eq. (|2"2"|) accounts for the response of the subsystem n to a force acting on the DoF next to the 
(np) interface. We define the part of the Green function G^ np ^ which gives the wave field on an (np 1 ) interface due to an 
excitation on the (np) interface as G^ pnp '. The sub-matrices G tyPnp ' can now be calculated by FEM. 



4.2. Flux formulas 

The energy flux from the (np) to the (np') interface produced by the force f in the n t h deterministic subsystem can be 
written in a form similar to Eq. (II9p . that is, 

P i%' = |s{(f (np) ) + (G^ np,) ) + ^ {d[7' ] } G ( d pnp ' ) i {np) }- (23) 

One can also find the total flux of energy flowing from the pth stochastic subsystem into the nth deterministic subsystem, 
that is, 

Pp ^ n = |3{{q™ + G (pnp) t [np) y^ {D% P } } G {pnp) t ("*>}. (24) 

Note that 

p p 
p—i p,p'—i 

in general; equality holds only if there is no energy loss due to, for example, damping in subsystem n. 

So far everything is exact; however, the stochastic force given in Eq. (|2ip is not known as in depends on qi„. We proceed 
now by taking averages over the power inputs in terms of the reverberant field components. That is, we consider averages (P) 
by determining fluxes over an ensemble of stochastic subsystems where we make small changes to the random boundaries. 
Using Eqs. (|2"Tj) and ([2"5]h one can write 

( p pt) - f s{E^ ( r p>) <^jv„>}, (25) 

where 

T (p^ P ',n) = D t. {GfV^sft {d&*">} G {pnp,) D dh (26) 
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where, according to Eq. (|24l) . 



and indexes j, k run over DoFs on the corresponding interface. The average energy influx from the pth to the nth subsystem 
can be written in the form 

= ^{J2 T lk n ((iin)*ii)h (27) 

(l + DUGt^f) ® {d ( 7 ] } G ( r P) D dl . (28) 

Here, the identity matrix accounts for the incoming wave energy, whereas outgoing flux on the (np) interface due to reflections 
inside the nth subsystems is described by the Green function G^ pnp ' . 

4.3. The field-field correlation function 

The field variables qj„ enter now only in form of the correlation functions {{q{ n ) Qi n }- The latter can be express in the 
following form 

((qiSqt)=C 2 p F(r J ,r k ), (29) 

where F(r 3 , r^) is the area (volume) normalized field amplitude correlation function, and C p is the mean squared amplitude 
of the incoming component of the random wavefield in the pth subsystem. Taking into account that incoming and outgoing 
components (with respect to any direction) of the diffuse field store the same amount of energy, one can find the time and 
ensemble averaged kinetic energy in the pth subsystem as 

(4 km) ) = (30) 

where p p is the material density and S p is the area (volume) of the subsystem. Finally, using the virial theorem, one can 
express C p through the total energy 

The problem is thus reduced to finding the correlation function F{Tj,Tk). The issue of random field correlation function 
was discussed in [15] in the context of random diffuse field theory or in [23 in the context of quantum and wave chaos. The 
central result is that the correlation function for a diffuse field is (again after averaging over an ensemble of similar systems) 
equal to the imaginary part of the Green function of this system [2~il [2"5]. In particular, the area- normalized correlation 
function is given in the following form 

F ( r ^ fc ^ S{ tr( G t,r fc ))} 5{G(r - rfc)} ' (32) 

where G(rj, r^) is he Green function connecting points i"j,rfc and trG denotes the trace. Taking into account that the trace 
of the Green function is related to the modal density, one can write 

„( w ) = -^9f{tr(G(r J -,r fc ))}. (33) 

The modal density of the subsystem is to leading order given by the area (volume) of the subsystem; in the two dimensional 
case the expression for the modal density is thus of the form 

"H - fg, (34) 

where c is the speed of sound. Corrections to the leading order behavior for the modal density due to, for example, the 
boundary of the subsystem can be considered in terms of a Weyl expansion, see [TBI [23] ■ Finally, using Eqs. ((32]) . (|33|) . one 
finds 

<(«DV n > = -^r^Gfe,^)}. (35) 
nujn(uj ) 

It may seem that the result obtained does not simplify the problem, because we do not know the Green function for the 
whole system. However, we can approximate the global solution locally using the free-space Green function. In fact, this 
issue was thoroughly investigated in wave chaos theory, see [TB, 26 for a reviews. It turns out that the imaginary part of the 
free-space Green function describes the correlations in Gaussian random fields correctly. The free-space Green function can 
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indeed be viewed as the zeroth-order approximation to the true correlation function |27j . For the 2D Hclmholtz equation 
considered here, the free Green function is given as 

G (|r|) = iy (7i|r|) - | J (7i|r|)}, (36) 

where the complex part of z = —irjeu + uj 2 takes into account absorption. In the case r\ = and using the free Green function, 
yields the well known result for Gaussian random waves |28j 

F(|r|) = JoMr|). 

The average flux in Eqs. (|25p and (|27[) is found by using a uniform diffuse field assumption with uniform reverberant field 
amplitudes C p across the boundary. This is a basic assumption of SEA and the coefficients T^ p ^ rp can be related to SEA 
coupling constants. Moreover, we assumed that the diffuse field at different interfaces of a given deterministic subsystem 
are not correlated. The random boundaries in each subsystem are indeed assumed to vary independently; if furthermore 
the lengths of the random boundaries of the stochastic subsystems exceed the lengths of the interfaces, the reverberant field 
in each subsystem is most strongly influenced by its boundaries. Thus, one can neglect coherent power transfer between 
two stochastic subsystems 29J and the average net flux carried by the diffuse field through any section of a deterministic 
subsystem is the sum of contributions given by Eqs. (|23|) . (j2"7) . 



5. Energy balance equations 



We have now all the ingredients to determine the total energy influx into the pth stochastic subsystem through its interfaces, 
that is, 

1 p A p AE {r , ] ) 

m=^+h^-fir-' (37) 

p'=l " p p 

where the direct field contribution Pp d ^ is given in (|19D and 

N 

^ = E S KE r i¥" tf,W)f, ('J» r *)}5 (38) 

n=l j,k 

note, that A^, is zero if the p' th subsystem is not directly connected to the p t h subsystem via a deterministic subsystem. 

The mean energy flowing out of the pth subsystem either due to damping and or due to energy transfer into neighboring 
deterministic subsystems can be written as 

Sp 

where 

^-x^E^r"^^)} (40) 

n=l j,k 

and r)p is the damping factor. We assume here that rj p is constant across each stochastic subsystem. The first term in 
Eq. (|3T))) is the energy lost due to damping in the direct field which can be found with the help of appropriate direct 
field propagators L. The second term is the energy lost in subsystem p due to damping and the last term gives the energy 
transferred into adjacent deterministic subsystems. 

We finally obtain a set of SEA-like equations using the energy balance relations (Pp n ) = (P° ut ). That makes it possible 
to reduce the problem to a set of linear equations 

(Pf n ) - (P? ut ) = 

(Pt) - (pr*> = o (41) 

(Pp n ) - (P£ ut ) = o 

which can be solved for the unknowns (E[ r) ), (E { 2 r) ),. . . , (E { p ] ). The energies stored in the direct field (E^), . . . , can 
again be found with the help of the direct field propagator L. 
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We now rewrite the set of equations in matrix form, that is, 

/ M n - m - 

pi 

\ M lP 

where the diagonal terms M pp are given as 



the off-diagonal terms M pp > have the form 



and the source terms are 

Q. 

We can now directly read off coupling loss factors and similar SEA parameters from Eq. (J22J) - Note that in contrast to 
traditional SEA, we can now incorporate losses due to damping within the deterministic subsystems and the summation over 
the row vectors in (|42[) may not add up to zero. 

One can also find the reverberant energies in the n t h deterministic subsystem as an incoherent sum of reverberant 
contribution from all neighboring subsystems 

p 

where 

< E S>> = ^S)S^' ,3{G ° <ri ' r * )} ' (47) 

and 

= Dj(G[ np) )^ [D^G^ p) D dl . (48) 

Our approach is thus a variant of SEA with respect to treating the stochastic components. The coupling loss factors 
contain wave information from the deterministic subsystems - they keep track of resonances excited in the deterministic 
subsystems, but coupled to the stochastic components. In particular, our hybrid method provides 

(i) direct coupling between deterministic subsystems via the direct field which is here calculated explicitly and coherently 
across the whole structure - a contribution not considered so far in hybrid approaches. 

(ii) an exact numerical treatment of the deterministic subsystems giving rise to coupling between stochastic subsystems. 
The method thus yields frequency dependent SEA coupling parameters. 

(iii) a statistical treatment of the stochastic subsystems. Using the correlation function 

F(T j ,T k )={p(T j )4,(T k )) 

is equivalent to the assumption that the underlying reverberant field is diffuse. Our approach is equivalent to using the direct 
field reciprocity relation derived in |13j : it gives, however, a more straightforward access to the problem. 

6. Numerical validation 

The method has been numerically validated for the structure shown in Fig. [T] with a 6- function force in one of the deterministic 
subsystem as shown in Fig. [31 We solve the wave problem using both the hybrid approach and by directly applying the 
FEM. To introduce randomness into the latter, we slightly distort the shape of the boundaries in the stochastic subsystems 
keeping the areas constant and find the FEM solution for every realization. Depending on the damping parameter, there are 
two main scenarios. In the high damping case, the energy flux into the stochastic subsystems is damped before it reaches 
the random boundaries and reflections into the reverberant field are suppressed. A large part of the energy is thus stored 
in the direct field. In the low damping case, energy loss in the direct field is insignificant compared to the net flux into the 
stochastic components and most of the energy is stored in the stochastic reverberant field. In this case, almost all waves 





(42) 



M pp = 



M pp , = 



1 Bp 

2,uj PpSp 

i K> 

2,lu Ppi Spi 



(43) 
(44) 
(45) 
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Figure 4. The energies in the left (pi) and right (^2) region found with FEM using Monte Carlo sampling over randomized 
boundaries compared to the results obtained with the hybrid method. Light grey: each Monte Carlo realization; black: 
ensemble average; bold grey: hybrid method; 77 = 0.20. 



reach the random boundaries and are reflected into the reverberant field where they accumulate. Increasing the damping 
parameter 77, one can cross over from one to the other regime. 

In Fig. we present results for the high damping case. The response was found as an average over ten realizations of the 
random boundaries. One can see a good agreement between the hybrid result and the full FEM calculation in both stochastic 
subsystems. Note the difference in scale for subsystem p\ and pi\ as pi is next to the deterministic subsystem containing 
the source and damping is relatively large, most of the energy is stored in this subsystem. In addition, more than 80% of 
the energy is stored in the direct field over the whole frequency range, see Fig. [5] The response of each individual realisation 
of the ensemble thus coincides well with the ensemble average as shown in Fig. UJ The signal in p\ , on the other hand, 
fluctuates wildly over the different ensemble realisations. The hybrid method yields the ensemble averaged response and the 
amount of energy in the direct and reverberant field components are comparable. In Fig. [5j we plot the direct field energies 
Epf 1 , Epf in regions p\ and p 2 together with the total energies (E Pl ) , (E P2 ) as function of the frequency. The resonance 
features of the response in Fig. [5] are mostly due to resonances in the subsystem ri2- Note that the resonance positions in 
both field components roughly coincide, although small shifts are observable. The relative intensity of resonance peaks is 
also different for the two field components. Clearly, the reverberant and direct field couple differently with the resonances in 
the deterministic subsystems. 

In the high damping regime, the frequency dependence of the response in the FEM calculations is relatively smooth 
because of the high modal overlap, see Fig. SJ That is, individual resonances in the stochastic subsystems tend to be broad 
and are not resolved. In contrast, in the low damping regime, one observes sharp peaks in each individual FEM realization 
which are associated with individual resonances of the stochastic subsystems. As a consequence, the response function 
varies considerably across the ensemble and individual realizations have little in common apart from the resonance structure 
originating from the deterministic subsystems. The computational results are presented in Fig. [5] where one can clearly see 
the response variation across the ensemble. (Note that results are shown on a logarithmic scale). The response was found 
as an average over 36 realizations. Again, one observes a good agreement between the hybrid method and the ensemble 
average. In Fig. [6l we also present the results for the mean energy obtained by conventional SEA. One can see that SEA is 
unable to resolve any features due to resonant properties of the deterministic subsystems. It just gives a ballpark figure for 
the expected wave intensities which is a constant for our damping model (with damping parameter 77 independent of u>) . 
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Figure 5. (Colour online) Light area: the mean total energy in the left (pi) and right ($2) plates; dark area: the energy in 
the direct field Epf* . The inset in the left figure shows the ratio of the energy in the direct field Ep\ ' to the total energy E P1 



7. Conclusion 



A hybrid FEM/SEA method has been developed which works efficiently in the mid- frequency regime. The key features of 
the method are an exact computation of the direct field over the whole system and the use of random diffuse field correlation 
functions to evaluate the coupling constants between reverberant field components corresponding to different stochastic 
subsystems. The direct and reverberant fields are then coupled through energy balance equations after FEM models for 
both the reverberant and the direct field have been solved. The FEM model for the direct field acts globally across the 
system coupling local deterministic models via direct field propagators. For the stochastic components, we solve the wave 
problem locally by finding Green functions for the deterministic subsystem. Our method has been tested in both the high 
and low damping regime. For high damping, a large part of the energy is stored and dissipated in the direct field. The 
energy distribution in the direct field is highly non-uniform within each stochastic subsystem and coupling between different 
stochastic subsystems can therefore not be predicted by SEA [30] . The method developed here makes it possible to overcome 
this difficulty and introduces coherent coupling mechanism between different deterministic subsystems. It is valid for all 
penetration depths of acoustic vibrations into the medium and is able to describe the cross-over from the high to the low 
damping regime when most of energy is stored in the reverberant field. It, furthermore, resolves resonance structures due to 
the deterministic subsystems in both the direct and reverberant field. 

Like the hybrid method developed in 112) , we determine local coupling parameters between stochastic subsystems 
using FEM which takes into account excitations through a diffuse field in the stochastic subsystems. In contrast to Shorter 
and Langley, our approach makes directly use of the result Eq. (l3"5j) for the field-field correlation function instead of using 
a reciprocity relationships between direct field radiation and diffuse reverberant loading |13) . In addition, we calculate 
the direct field explicitly which is particularly important in the high damping limit. Thus, the method is applicable for 
over-damped structures where pure SEA may not produce a reliable result (31 j - 
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Figure 6. The low damping regime: energies in the left and right region found with Monte Carlo sampling compared to 
the results obtained with the hybrid method. Light grey: each Monte Carlo realization; black: ensemble average; bold grey: 
hybrid method; dashed line the SEA result; r\ — 0.01. 
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